MSIS Version Evaluation- Multiple Arcs¶
Details:¶
This notebook shows GEODYN II output for from three versions of the MSIS models for thermospheric density: MSISe86, MSISe00, MSISe2.0. The beginning of the notebook contains the analysis that can be gathered from assessing the plots that are demonstrated below.
Run Parameters:
Satellite: Starlette Observation Datatype: SLR Accel9 Adjustments: Off Density Models: NRLMSISe86, NRLMSISe00, MSIS 2.0
Data arcs in this analysis:
Arc 1 -------- 14 days ------ 03/09/14 - 03/09/28 Arc 2 -------- 14 days ------ 03/09/28 - 03/10/12 Arc 3 -------- 14 days ------ 03/10/12 - 03/10/26 Arc 4 -------- 14 days ------ 03/10/26 - 03/11/09 Arc 5 -------- 14 days ------ 03/11/09 - 03/11/23 Arc 6 -------- 14 days ------ 03/11/23 - 03/12/07 Arc 7 -------- 14 days ------ 03/12/07 - 03/12/21 Arc 8 -------- 14 days ------ 03/12/21 - 04/01/5
## Analysis |
### Background: |
Starlette Orbit Info: |
|
Increased density from Geomagnetic and solar activity : |
the increase in solar flux than due to increased Ap. |
### Relative model performance: |
We have to look at a combination of different factors. This notebook includes the following: 1. RMS of fit on final iteration (want small) 2. Cd corrections from a priori (want small) - In general, the drag coefficient will be adjusted (corrected by some amount) to account for variations seen in the mass density provided by the different models. This makes the Cd (in this run scenario) a physically meaningful parameter to assess the POD. 3. Observation residuals (want to be small, but don’t rely too heavily on this) |
The mean in the residuals (or the mean by station) is not important, it is very small in these runs |
Adjusted Cd (apriori value = 2.2) - At this altitude regime (and over this epoch), the MSISe86 model has the smallest corrections from a priori at a mean percent difference of 18.4%. MSISe00 has slightly higher corrections at 19.1 %. MSISe2.0 has a mean percent difference from a priori of 20.25%. - In stormtime scenarios, the Cd respond by correcting the mean density provided by the different models. - At this altitude regime, the MSIS2.0 has the largest correction during storm time events. The end of Arc 5 shows a particularly large correction in which all models correct the Cd by more than 200% (m2=~280%, m00=~240%, m86=~205%). - from a POD perspective, the Cd does a good job of accounting for the density variations. Its not clear to me at this stage what the physical significance of the corrections is… [ASK QUESTION] |
RMS of Fit |
|
It is worth noting that since this is at such a high altitude, the many improvements that were made to the MSIS models in version 00 and 2.0 may not be impacting these results in a POD sense. The data that went into MSIS2.0 was mostly from satellites and data sets derived at lower altitudes. We will need to repeat the assessment methods demonstrated here for a lower altitude satellite geodyn run and compare the results. |
Load the simulations over all arcs¶
[1]:
# Import the geodyn output readers
import sys
sys.path.insert(0, '/data/geodyn_proj/analysis/util_funcs/py_read_geodyn_output/')
from a_ReadGEODYN import Read_GEODYN_func_MultipleArcs
[2]:
%load_ext autoreload
%autoreload 2
# Load multiple arcs:
arc_list = ['030914_2wk',
'030928_2wk',
'031012_2wk',
'031026_2wk',
'031109_2wk',
'031123_2wk',
'031207_2wk',
# '031221_2wk', # this arc has some errors at the moment...
]
geodyn_run_params_m1 = {}
geodyn_run_params_m1['arc_list'] = arc_list
geodyn_run_params_m1['sat_file'] = 'st'
geodyn_run_params_m1['SAT_ID'] = 7501001
geodyn_run_params_m1['grav_id'] = 'goco05s'
geodyn_run_params_m1['local_path'] = '/data/geodyn_proj/analysis/starlette_analysis/'
geodyn_run_params_m1['AccelStatus'] = 'acceloff'
geodyn_run_params_m1['den_model'] = 'msis86'
path_to_model = ('/data/geodyn_proj/runs_geodyn/results/'+
geodyn_run_params_m1['sat_file'] +'/'+
geodyn_run_params_m1['den_model']+'/'+
geodyn_run_params_m1['den_model']+'_'+
geodyn_run_params_m1['AccelStatus'] +'/')
geodyn_run_params_m1['path_to_data'] = path_to_model
geodyn_run_params_m1['YR'] = 2003
geodyn_run_params_m1['DATA_TYPE'] = 'SLR'
geodyn_run_params_m1['verbose_loading'] = False
geodyn_run_params_m1['Verbose_Stats'] = False
[3]:
(AdjustedParams_m1,
trajectory_df_m1,
Density_m1,
resids_observed_m1,
resids_summary_station_m1,
resids_meas_summary_m1,
RunStats_m1) = Read_GEODYN_func_MultipleArcs(geodyn_run_params_m1)
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st030914_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st030928_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st031012_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st031026_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st031109_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st031123_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis86/msis86_acceloff/st031207_2wk.goco05s
[4]:
geodyn_run_params_m2 = geodyn_run_params_m1
geodyn_run_params_m2['den_model'] = 'msis00'
path_to_model = ('/data/geodyn_proj/runs_geodyn/results/'+
geodyn_run_params_m2['sat_file'] +'/'+
geodyn_run_params_m2['den_model']+'/'+
geodyn_run_params_m2['den_model']+'_'+
geodyn_run_params_m2['AccelStatus'] +'/')
geodyn_run_params_m2['path_to_data'] = path_to_model
(AdjustedParams_m2,
trajectory_df_m2,
Density_m2,
resids_observed_m2,
resids_summary_station_m2,
resids_meas_summary_m2,
RunStats_m2) = Read_GEODYN_func_MultipleArcs(geodyn_run_params_m2)
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st030914_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st030928_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st031012_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st031026_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st031109_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st031123_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis00/msis00_acceloff/st031207_2wk.goco05s
[5]:
geodyn_run_params_m3 = geodyn_run_params_m1
geodyn_run_params_m3['den_model'] = 'msis2'
path_to_model = ('/data/geodyn_proj/runs_geodyn/results/'+
geodyn_run_params_m3['sat_file'] +'/'+
geodyn_run_params_m3['den_model']+'/'+
geodyn_run_params_m3['den_model']+'_'+
geodyn_run_params_m3['AccelStatus'] +'/')
geodyn_run_params_m3['path_to_data'] = path_to_model
(AdjustedParams_m3,
trajectory_df_m3,
Density_m3,
resids_observed_m3,
resids_summary_station_m3,
resids_meas_summary_m3,
RunStats_m3) = Read_GEODYN_func_MultipleArcs(geodyn_run_params_m3)
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st030914_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st030928_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st031012_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st031026_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st031109_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st031123_2wk.goco05s
Loading /data/geodyn_proj/runs_geodyn/results/st/msis2/msis2_acceloff/st031207_2wk.goco05s
Import some helpful plotting schemes
[6]:
%load_ext autoreload
%autoreload 2
sys.path.insert(0, '/data/geodyn_proj/analysis/util_funcs/util_common/')
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
from plotly.subplots import make_subplots
import plotly.express as px
from plotting_scheme_functions import help_plot_colors, add_arc_background_w_text, legend_as_annotation
col1, col2, col3 = help_plot_colors()
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Coefficient of Drag Adjustments¶
Details and Background:¶
the \(C_d\) is an adjusted parameter in GEODYN. Each iteration, it is estimated and adjusted to improve the OD simulation
\(C_d\) is treated as a constant in GEODYN (for any iteration). The factor \(C_d\) varies slightly with satellite shape and atmospheric composition. However, for any geodetically useful satellite, it may be treated as a satellite dependent constant. If the Cd is time dependent, then it is constant over the designated time intervals.
$ \bar`{A}_D = -:nbsphinx-math:frac{1}{2}` C_d \frac{A_s}{m_s} \rho_d v_r :nbsphinx-math:`bar`{v}_r $
\(A_D\) acceleration due to atmospheric drag force
\(C_d\) is the satellite drag coefficient
\(A_s\) is the cross-sectional area of the satellite
\(m_s\) is the mass of the satellite
\(\rho_d\)is the-density of the atmosphere at the satellite position
\(\bar{v}_r\) is the velocity vector of the satellite relative to the atmosphere, and
\(v_r\) is the magnitude of the velocity vector, vr .
Cd Analysis¶
A few key features can be seen in the above plot: 1. When there is some sort of storm event early in the arc, the Cd responds by correcting the mean density provided by the different models. 2. The MSISe00 corrections are less at this set of altitudes as opposed to those of m86 and m2.0. 3. Overall, though from a POD perspective, they seem to do a good job accounting for density variations. - Frank states that the onus is to show there is some real storm event going on over these several days.
[8]:
SAT_ID = geodyn_run_params_m1['SAT_ID']
which_stat = 'CURRENT_VALUE' # 2 is the current val
all_cd_m1 = []
all_dates_m1 = []
for ii,arc in enumerate(arc_list):
i_arc = ii+1
last_iter = list(AdjustedParams_m1[arc].keys())[-1]
labels = list(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys())
for i,val in enumerate(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys()):
all_cd_m1.append(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'][val][which_stat])
all_dates_m1.append(labels[i])
all_cd_m2 = []
all_dates_m2 = []
for ii,arc in enumerate(arc_list):
i_arc = ii+1
last_iter = list(AdjustedParams_m2[arc].keys())[-1]
labels = list(AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'].keys())
for i,val in enumerate(AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'].keys()):
all_cd_m2.append(AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'][val][which_stat])
all_dates_m2.append(labels[i])
all_cd_m3 = []
all_dates_m3 = []
for ii,arc in enumerate(arc_list):
i_arc = ii+1
last_iter = list(AdjustedParams_m3[arc].keys())[-1]
labels = list(AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'].keys())
for i,val in enumerate(AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'].keys()):
all_cd_m3.append(AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'][val][which_stat])
all_dates_m3.append(labels[i])
# take % difference from a priori
cd_apriori= 2.2
percdiff_cd_m1 = (np.abs(np.array(all_cd_m1) - cd_apriori)/ cd_apriori)*100
percdiff_cd_m2 = (np.abs(np.array(all_cd_m2) - cd_apriori)/ cd_apriori)*100
percdiff_cd_m3 = (np.abs(np.array(all_cd_m3) - cd_apriori)/ cd_apriori)*100
avg_percdiff_cd_m1 = "{:.5g}".format(np.mean(percdiff_cd_m1))
avg_percdiff_cd_m2 = "{:.5g}".format(np.mean(percdiff_cd_m2))
avg_percdiff_cd_m3 = "{:.5g}".format(np.mean(percdiff_cd_m3))
print('Avg % diff for MSISe86: ',avg_percdiff_cd_m1,'%')
print('Avg % diff for MSISe00: ',avg_percdiff_cd_m2,'%')
print('Avg % diff for MSISe2.0:',avg_percdiff_cd_m3,'%')
Avg % diff for MSISe86: 18.466 %
Avg % diff for MSISe00: 19.128 %
Avg % diff for MSISe2.0: 20.251 %
Adjusted Cd Plot (timeseries)¶
[9]:
SAT_ID = geodyn_run_params_m1['SAT_ID']
which_stat = 'CURRENT_VALUE' # 2 is the current val
# make plot:
fig = make_subplots(
rows=2, cols=1,
subplot_titles=(["Timeseries of Cd", "Percent difference from a priori (Cd=2.2)"]),
vertical_spacing = 0.1,
)
for ii,arc in enumerate(geodyn_run_params_m1['arc_list']):
i_arc = ii+1
last_iter = list(AdjustedParams_m1[arc].keys())[-1]
labels = list(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys())
val_list_1 = []
for i in AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_1.append(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_1,
name= geodyn_run_params_m1['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[0],
size=5,),
showlegend=False,
),
row=1, col=1,
)
last_iter = list(AdjustedParams_m2[arc].keys())[-1]
val_list_2 = []
for i in AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_2.append(AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_2,
name=geodyn_run_params_m2['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[1],
symbol = 'circle',
size=5,),
showlegend=False,
),
row=1, col=1,
)
last_iter = list(AdjustedParams_m3[arc].keys())[-1]
val_list_3 = []
for i in AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_3.append(AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_3,
name= geodyn_run_params_m3['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[2],
symbol = 'circle',
size=5,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=all_dates_m1,
y=percdiff_cd_m1,
name= 'MSIS 86',
mode='markers',
marker=dict(
color=col1,
size=5,
),
showlegend=False,
),
row=2, col=1,)
fig.add_trace(go.Scatter(x=all_dates_m2,
y=percdiff_cd_m2,
name= 'MSIS 00',
mode='markers',
marker=dict(
color=col2,
size=5,
),
showlegend=False,
),
row=2, col=1,)
fig.add_trace(go.Scatter(x=all_dates_m3,
y=percdiff_cd_m3,
name= 'MSIS 2',
mode='markers',
marker=dict(
color=col3,
size=5,
),
showlegend=False,
),
row=2, col=1,)
fig.add_annotation(
x=.95,
y=.98,
xref="x domain",
yref="y domain",
showarrow=False,
text='msis 86: '+avg_percdiff_cd_m1 +'%',
font=dict(
size=14,
color="#ffffff"
),
align="center",
bordercolor="#c7c7c7",
borderwidth=2,
borderpad=4,
bgcolor=col1,
opacity=0.9,
row=2, col=1,
)
fig.add_annotation(
x=.95,
y=.88,
xref="x domain",
yref="y domain",
showarrow=False,
text='msis 00: '+avg_percdiff_cd_m2 +'%',
font=dict(
size=14,
color="#ffffff"
),
align="center",
bordercolor="#c7c7c7",
borderwidth=2,
borderpad=4,
bgcolor=col2,
opacity=0.9,
row=2, col=1,
)
fig.add_annotation(
x=.95,
y=.78,
xref="x domain",
yref="y domain",
showarrow=False,
text='msis 2: '+avg_percdiff_cd_m3 +'%',
font=dict(
size=14,
color="#ffffff"
),
align="center",
bordercolor="#c7c7c7",
borderwidth=2,
borderpad=4,
bgcolor=col3,
opacity=0.9,
row=2, col=1,
)
fig = add_arc_background_w_text(fig, 9.5, arc_list, True)
fig = legend_as_annotation(fig)
# fix layout info:
fig.update_yaxes( title="Cd ",exponentformat= 'power',row=1, col=1)
fig.update_yaxes( title="% difference",exponentformat= 'power',row=2, col=1)
fig.update_xaxes( title="Date", row=2, col=1)
fig.update_layout(title="Time Dependent Drag Coefficient ")
fig.update_layout(
font=dict(
size=14,
),)
fig.update_layout(
autosize=False,
width=900,
height=900,
)
iplot(fig)
Density Output¶
The density output is sourced from the
DRAG.f90subroutine.
Density output Analysis¶
Starlette is in somewhat of an eccentric orbit, roughly 800 km x 1100 km orbit - so that accounts for the change in density along the orbit.
while this is helpful to see, one really needs to demonstrate these results over multiple arcs.
The relative densities are different at Apogee and perigee.
Apogee : MSISe86 and MSISe00 have the higher densities
Perigee : MSIS2.0 has the higher density
[10]:
# make some space by deleting the columns that wont be used here
del_columnns = ['year',
'month',
'day',
'hours',
'minutes',
'secs',
'timeHHMMSS',
'Elapsed',
'YYMMDD',
'HHMMSS']
for arc in arc_list:
for i in del_columnns:
if i in Density_m1[arc]:
del Density_m1[arc][i]
if i in Density_m2[arc]:
del Density_m2[arc][i]
if i in Density_m3[arc]:
del Density_m3[arc][i]
Density timeseries plot¶
[11]:
data_nums_1 = 1000
data_nums_2 = 200
fig = make_subplots(rows=1, cols=1,
subplot_titles=(["Every 200th point"]),
)
for arc in arc_list:
fig.add_trace(go.Scatter(x=Density_m1[arc]['Date'][::data_nums_2],
y=Density_m1[arc]['rho (kg/m**3)'][::data_nums_2],
name = geodyn_run_params_m1['den_model'],
mode='markers',
marker=dict(color=col1,
size=3,
),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=Density_m2[arc]['Date'][::data_nums_2],
y=Density_m2[arc]['rho (kg/m**3)'][::data_nums_2],
name = geodyn_run_params_m2['den_model'],
mode='markers',
marker=dict(color=col2,
size=3,
),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=Density_m3[arc]['Date'][::data_nums_2],
y=Density_m3[arc]['rho (kg/m**3)'][::data_nums_2],
name = geodyn_run_params_m3['den_model'],
mode='markers',
marker=dict(color=col3,
size=3,
),
showlegend=False,
),
row=1, col=1,
)
fig.update_layout(
title="Density along Starlette Orbit",
)
fig.update_yaxes(type="log", exponentformat= 'power',)
fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_yaxes(title_text="kg/m^3", row=1, col=1)
fig.update_xaxes(title_text="Date", row=1, col=1)
fig.update_layout(legend= {'itemsizing': 'constant'})
# fig = add_arc_background_w_text(9.8e-14)
# fig = legend_as_annotation()
fig = add_arc_background_w_text(fig, 9.8e-14, arc_list, True)
fig = legend_as_annotation(fig)
fig.update_layout(
autosize=False,
width=900,
height=700,
)
fig.update_layout(
font=dict(
size=18,
),)
iplot(fig)
[12]:
# load the saved F107 and AP for these arcs from a saved file.
f107a =[]
f107d =[]
date = []
ap = []
import pandas as pd
for i,arc in enumerate(arc_list):
# print(i,arc)
msisin_den_file = geodyn_run_params_m3['path_to_data'] + "DENSITY/st"+arc + ".goco05s_msisin"
DEN1_csv = pd.read_csv(msisin_den_file,
skiprows = 1,
names = ['IYYDDD',
'IYR',
'DAY',
'UTSEC',
'ALTKM',
'GLAT',
'GLON',
'STLOC',
'AVGFLX',
'FLUX',
'AP1',
'AP2',
'AP3',
'AP4',
'AP5',
'AP6',
'AP7',
],
sep = '\s+',
)
DEN1_csv['Date'] = (pd.to_datetime('0'+ ((DEN1_csv['IYR'].astype(int).astype(str))), format='%y')
+ pd.to_timedelta(DEN1_csv['DAY'], unit='days'))
f107a.append(DEN1_csv['AVGFLX'].values)
f107d.append(DEN1_csv['FLUX'].values)
date.append(DEN1_csv['Date'].values)
ap.append(DEN1_csv['AP7'].values)
[13]:
# import a function that will construct the orbit averaged density
from orbit_average_func import orb_avg
Orbit averaged Density Plot¶
[14]:
fig = make_subplots(rows=3, cols=1,
subplot_titles=(['Orbit Averaged Density', 'Ap', 'F10.7']),
horizontal_spacing = 0.05,
# vertical_spacing = 0.12,
vertical_spacing = 0.1,
)
for i_arc,arc in enumerate(arc_list):
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m1, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col1,
size=3,),
showlegend=False,
),
row=1, col=1,
)
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m2, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col2,
size=3,),
showlegend=False,
),
row=1, col=1,
)
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m3, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col3,
size=3,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=pd.to_datetime(date[i_arc][::500]),
y=ap[i_arc][::500],
name= 'Arc :'+str(i_arc+1),
mode='markers',
marker=dict(color = 'Grey',
size=3,),
showlegend=False,
),
row=2, col=1,
)
fig.add_trace(go.Scatter(x=pd.to_datetime(date[i_arc][::1000]),
y=f107d[i_arc][::1000],
# name= 'Ap',
mode='markers',
marker=dict(color = 'Grey',
size=3,),
showlegend=False,
),
row=3, col=1,
)
fig = add_arc_background_w_text(fig, 1.1*np.max(d_avg_rolling ), arc_list, False)
fig = legend_as_annotation(fig)
fig.update_layout(
title="Density along Starlette Orbit",
)
fig.update_yaxes(type="log", exponentformat= 'power',row=1, col=1)
# fig.update_xaxes(title_text="Date", row=1, col=1)
# fig.update_xaxes(title_text="Date", row=2, col=1)
fig.update_xaxes(title_text="Date", row=3, col=1)
fig.update_yaxes(title_text="kg/m^3", row=1, col=1)
fig.update_yaxes(title_text="nT", row=2, col=1)
fig.update_yaxes(title_text="sfu", row=3, col=1)
fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_layout(
autosize=False,
width=900,
height=1000,)
fig.update_layout(
font=dict(
size=18,
),)
iplot(fig)
Composite Plots, density, cd, activity indicies¶
[15]:
fig = make_subplots(rows=3, cols=1,
subplot_titles=(['Orbit Averaged Density','Adjusted Cd', 'Geomagnetic and Solar Activity']),
horizontal_spacing = 0.05,
# vertical_spacing = 0.12,
vertical_spacing = 0.1,
specs=[[{"secondary_y": False}],
[{"secondary_y": False}],
[{"secondary_y": True}], ])
for i_arc,arc in enumerate(arc_list):
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m1, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col1,
size=3,),
showlegend=False,
),
row=1, col=1,
)
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m2, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col2,
size=3,),
showlegend=False,
),
row=1, col=1,
)
time_avg,d_avg, d_avg_rolling = orb_avg(Density_m3, arc)
fig.add_trace(go.Scatter(x=time_avg[::2],
y=d_avg_rolling,
name= ' Arc ' +str(i_arc) ,
mode='markers',
marker=dict(
color=col3,
size=3,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=pd.to_datetime(date[i_arc][::500]),
y=ap[i_arc][::500],
name= 'Ap',
mode='markers',
marker=dict(color = 'blue',
size=3,),
showlegend=False,
),
secondary_y=False,
row=3, col=1,
)
fig.add_trace(go.Scatter(x=pd.to_datetime(date[i_arc][::1000]),
y=f107d[i_arc][::1000],
name= 'F10.7',
mode='markers',
marker=dict(color = 'red',
size=3,),
showlegend=False,
),
secondary_y=True,
row=3, col=1,
)
SAT_ID = geodyn_run_params_m1['SAT_ID']
which_stat = 'CURRENT_VALUE' # 2 is the current val
for ii,arc in enumerate(geodyn_run_params_m1['arc_list']):
i_arc = ii+1
last_iter = list(AdjustedParams_m1[arc].keys())[-1]
labels = list(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys())
val_list_1 = []
for i in AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_1.append(AdjustedParams_m1[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_1,
name= geodyn_run_params_m1['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[0],
size=7,),
showlegend=False,
),
row=2, col=1,
)
last_iter = list(AdjustedParams_m2[arc].keys())[-1]
val_list_2 = []
for i in AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_2.append(AdjustedParams_m2[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_2,
name=geodyn_run_params_m2['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[1],
symbol = 'circle',
size=7,),
showlegend=False,
),
row=2, col=1,
)
last_iter = list(AdjustedParams_m3[arc].keys())[-1]
val_list_3 = []
for i in AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'].keys():
val_list_3.append(AdjustedParams_m3[arc][last_iter][SAT_ID]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
y=val_list_3,
name= geodyn_run_params_m3['den_model'] + ' | Arc ' +str(i_arc) +' | Iters: '+ str(last_iter),
mode='markers',
marker=dict(
color=px.colors.qualitative.Plotly[2],
symbol = 'circle',
size=7,),
showlegend=False,
),
row=2, col=1,
)
# fix layout info:
fig.update_yaxes( title="Cd ",exponentformat= 'power',row=2, col=1)
fig = legend_as_annotation(fig)
fig = add_arc_background_w_text(fig, 1.1*np.max(d_avg_rolling ), arc_list, False)
fig.update_layout(
title="Density, Cd, and Activity Indicies",
)
fig.update_yaxes(type="log", exponentformat= 'power',row=1, col=1)
# fig.update_xaxes(title_text="Date", row=1, col=1)
# fig.update_xaxes(title_text="Date", row=2, col=1)
fig.update_xaxes(title_text="Date", row=3, col=1)
fig.update_yaxes(title_text="kg/m^3", row=1, col=1)
fig.update_yaxes(title_text="Ap [nT]", row=3, col=1, secondary_y=False, color='blue')
fig.update_yaxes(title_text="F10.7 [sfu]", row=3, col=1, secondary_y=True, color='red')
fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_layout(
autosize=False,
width=850,
height=950,)
fig.update_layout(
font=dict(
size=18,
),)
iplot(fig)
Residuals – obervations¶
Details:¶
the residuals being output by geodyn should be rescaled to CM instead of METER. Most of your stations have cm level (or better) data precision.
The mean in the residuals (or the mean by station) is not important - it is very small in these runs - unless you turn out to have a station with a large systematic range bias - due to some temporary anomaly.
If you look at the setups, you’ll see sets of “MBIAS” cards – these adjust different sets of range biases by station per arc or by station per data pass based on recommendations from the ILRS. They are adjusted as a normal part of the POD process.
Analysis¶
You see the bowtie effect in the residuals above. That is a standard effect of least squares - where fits seem the best in the center of the arc.
Observed Residuals timeseries plot¶
[16]:
fig = go.Figure()
for i,arc in enumerate(arc_list):
i_arc = i+1
fig.add_trace(go.Scatter(x=resids_observed_m1[arc]['Date'][::3],
y=resids_observed_m1[arc]['Residual'][::3].values.astype(float),
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col1,
size=3,),
showlegend=False,
),
)
fig.add_trace(go.Scatter(x=resids_observed_m2[arc]['Date'][::3],
y=resids_observed_m2[arc]['Residual'][::3].values.astype(float),
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col2,
size=3,),
showlegend=False,
),
)
fig.add_trace(go.Scatter(x=resids_observed_m3[arc]['Date'][::3],
y=resids_observed_m3[arc]['Residual'][::3].values.astype(float),
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col3,
size=3,),
showlegend=False,
),
)
fig = add_arc_background_w_text(fig, .7, arc_list, True)
fig = legend_as_annotation(fig)
fig.update_layout(
title='Observation Residuals',
yaxis_title="Residuals",
xaxis_title="Date",
)
fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_layout(
font=dict(
size=17,
),)
iplot(fig)
Residual Measurement Summary plots (RMS of fit)¶
[17]:
import datetime
mark_size = 10
fig = make_subplots(rows=2, cols=1,
subplot_titles=(["Mean Residuals", 'RMS of Fit']),
vertical_spacing = 0.1,
)
for i,arc in enumerate(geodyn_run_params_m1['arc_list']):
i_arc = i+1
arc_val = arc[:6]
arc_date = np.datetime64(datetime.datetime.strptime(arc_val, '%y%m%d'))
iter_index = int(resids_meas_summary_m1[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m1[arc]['MEAN'].values.astype(float)[iter_index]*1e2
rms_arc = resids_meas_summary_m1[arc]['RMS'].values.astype(float)[iter_index]
fig.add_trace(go.Scatter(x=[arc_date],
y=[mean_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col1,
size=10,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=[arc_date],
y=[rms_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col1,
size=10,),
showlegend=False,
),
row=2, col=1,
)
for i,arc in enumerate(geodyn_run_params_m2['arc_list']):
i_arc = i+1
arc_val = arc[:6]
arc_date = np.datetime64(datetime.datetime.strptime(arc_val, '%y%m%d'))
iter_index = int(resids_meas_summary_m2[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m2[arc]['MEAN'].values.astype(float)[iter_index]*1e2
rms_arc = resids_meas_summary_m2[arc]['RMS'].values.astype(float)[iter_index]
fig.add_trace(go.Scatter(x=[arc_date],
y=[mean_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col2,
size=9,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=[arc_date],
y=[rms_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col2,
size=9,),
showlegend=False,
),
row=2, col=1,
)
for i,arc in enumerate(geodyn_run_params_m3['arc_list']):
i_arc = i+1
arc_val = arc[:6]
arc_date = np.datetime64(datetime.datetime.strptime(arc_val, '%y%m%d'))
iter_index = int(resids_meas_summary_m3[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m3[arc]['MEAN'].values.astype(float)[iter_index]*1e2
rms_arc = resids_meas_summary_m3[arc]['RMS'].values.astype(float)[iter_index]
fig.add_trace(go.Scatter(x=[arc_date],
y=[mean_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col3,
size=8,),
showlegend=False,
),
row=1, col=1,
)
fig.add_trace(go.Scatter(x=[arc_date],
y=[rms_arc],
name= 'Arc: '+str(i_arc),
mode='markers',
marker=dict(color=col3,
size=8,),
showlegend=False,
),
row=2, col=1,
)
fig = add_arc_background_w_text(fig, .002, arc_list, False)
fig = legend_as_annotation(fig)
fig.update_layout(
title="Residual Measurement Summary",
)
fig.update_yaxes(title_text="Residual [cm]", row=1, col=1)
fig.update_yaxes(title_text="RMS", row=2, col=1)
fig.update_xaxes(title_text="Arc Dates", row=2, col=1)
# fig.update_xaxes(title_text="Arc", row=2, col=1)
fig.update_layout(legend= {'itemsizing': 'constant'})
fig.update_xaxes(tickangle=0)
fig.update_layout(
font=dict(
size=15,
),)
fig.update_layout(
autosize=False,
width= 750,
height= 750,
)
fig.show()
avg_rms_m1 = []
avg_resid_m1 = []
for i,arc in enumerate(geodyn_run_params_m1['arc_list']):
iter_index = int(resids_meas_summary_m1[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m1[arc]['MEAN'].values.astype(float)[iter_index]
rms_arc = resids_meas_summary_m1[arc]['RMS'].values.astype(float)[iter_index]
avg_rms_m1.append(rms_arc)
avg_resid_m1.append(mean_arc)
avg_rms_m2 = []
avg_resid_m2 = []
for i,arc in enumerate(geodyn_run_params_m2['arc_list']):
iter_index = int(resids_meas_summary_m2[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m2[arc]['MEAN'].values.astype(float)[iter_index]
rms_arc = resids_meas_summary_m2[arc]['RMS'].values.astype(float)[iter_index]
avg_rms_m2.append(rms_arc)
avg_resid_m2.append(mean_arc)
avg_rms_m3 = []
avg_resid_m3 = []
for i,arc in enumerate(geodyn_run_params_m3['arc_list']):
iter_index = int(resids_meas_summary_m3[arc]['Iter'].values[-1])-1
mean_arc = resids_meas_summary_m3[arc]['MEAN'].values.astype(float)[iter_index]
rms_arc = resids_meas_summary_m3[arc]['RMS'].values.astype(float)[iter_index]
avg_rms_m3.append(rms_arc)
avg_resid_m3.append(mean_arc)
avg_rms_m1 = "{:.5e}".format(np.mean(avg_rms_m1))
avg_resid_m1 = "{:.5e}".format(np.mean(avg_resid_m1)*1e2)
avg_rms_m2 = "{:.5e}".format(np.mean(avg_rms_m2) )
avg_resid_m2 = "{:.5e}".format(np.mean(avg_resid_m2)*1e2)
avg_rms_m3 = "{:.5e}".format(np.mean(avg_rms_m3) )
avg_resid_m3 = "{:.5e}".format(np.mean(avg_resid_m3)*1e2)
print('MSISe86 RMS ' ,avg_rms_m1)
print('MSISe00 RMS ' ,avg_rms_m2)
print('MSISe2 RMS ' ,avg_rms_m3)
print()
print('MSISe86 Resid ' ,avg_resid_m1)
print('MSISe00 Resid ' ,avg_resid_m2)
print('MSISe2 Resid ' ,avg_resid_m3)
MSISe86 RMS 9.22857e-02
MSISe00 RMS 9.23714e-02
MSISe2 RMS 9.31429e-02
MSISe86 Resid -3.28571e-02
MSISe00 Resid -3.28571e-02
MSISe2 Resid -3.71429e-02
Residuals – station summary¶
I have been told that these are not helpful for assessing density models.
unless to identify a really bad station (they can crop up from time to time because of a station anomaly). If one finds a bad station, then one can just delete the data. In this example, it looks like Riyadh (RIYA7832) seems to have a mean bias of 2 cm?. That is a bit large and might be a sign of station coordinate error or some other station-dependent error. One could solve for an MBIAS for that station over the arc to take care of that small problem - you would have to check that an MBIAS for that station is not already adjusted.
TODO: NEED TO ASK FL WHAT EXACTLY THESE DO
[ ]: